Data for testing IHS

library(tidyverse)
library(terra)
library(sf)
library(leaflet)
library(whitebox)
library(scico)

Calculating covariates

May as well start with the same set as the source papers - slope, curvatures, TWI ¯\_(ツ)_/¯

Using whitebox here as it’s geomorphometry algorithms apparently handle outliers/noise better in DEMs with projected coordinate systems vs the algorithms implemented in terra::terrain().

elev_dir <- file.path('data_spatial', 'elevation')
if(!dir.exists(file.path('data_spatial', 'terrain_morphometry'))) {
  dir.create(file.path('data_spatial', 'terrain_morphometry'))
}
terr_dir <- file.path('data_spatial', 'terrain_morphometry')

wbt_slope(
  dem    = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
  output = file.path(terr_dir, 'Slope_1m.tif'),
  units  = "degrees",
  wd     = getwd(),
  compress_rasters = TRUE)

wbt_plan_curvature(
  dem    = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
  output = file.path(terr_dir, 'PlanC_1m.tif'),
  log    = FALSE,
  wd     = getwd(),
  compress_rasters = TRUE)

wbt_profile_curvature(
  dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
  output = file.path(terr_dir, 'ProfC_1m.tif'),
  log = FALSE,
  wd = getwd(),
  compress_rasters = TRUE
)

slope <- rast(file.path(terr_dir, 'Slope_1m.tif'))
planc <- rast(file.path(terr_dir, 'PlanC_1m.tif'))
profc <- rast(file.path(terr_dir, 'ProfC_1m.tif'))

TWI is a little more complicated as the DEM needs to be hydrologically corrected first.

# breach depressions (less interferey than filling sinks)
wbt_breach_depressions_least_cost(
  dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
  output = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
  dist = 250,
  max_cost = NULL,
  min_dist = TRUE,
  flat_increment = NULL,
  fill = TRUE,
  wd = getwd(),
  compress_rasters = FALSE
)

dem <- rast(file.path(elev_dir, 'Lidar_FPS_1m.tif'))
dem_bd <- rast(file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'))

# check diff
dem_diff <- dem_bd - dem
dem_diff[dem_diff == 0] <- NA

Differences after breaching range between -2.21 m and 0.2 m and are largely confined to single-pixel channels in open water, so that’s fine.

# for TWI, slope needs to also come from the breached DEM
# for consistency
wbt_slope(
  dem    = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
  output = file.path(terr_dir, 'Slope_BD_1m.tif'),
  units  = "degrees",
  wd     = getwd(),
  compress_rasters = TRUE)

wbt_d_inf_flow_accumulation(
  input     = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
  output    = file.path(terr_dir, 'Dinf_SCA_1m.tif'),
  out_type  = "sca",
  threshold = NULL,
  log       = FALSE,
  clip      = FALSE,
  pntr      = FALSE,
  wd        = getwd(),
  compress_rasters = TRUE)

wbt_wetness_index(
  sca   = file.path(terr_dir, 'Dinf_SCA_1m.tif'),
  slope = file.path(terr_dir, 'Slope_BD_1m.tif'),
  output = file.path(terr_dir, 'TWI_1m.tif'),
  wd = getwd(),
  compress_rasters = TRUE
)

twi <- rast(file.path(terr_dir, 'TWI_1m.tif'))

Not so sure I like this output, but it gives us something to cluster…